#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <string.h> // memcpy
#include <fstream>
/*
  use M_PI const: 3.14159 
*/

// отрезок моделирования по X
#define A 0
#define B 2
#define SIGMA 1
#define K 1

// Определение количества шагов
int getStepsT(double deltaT, double time){
  return (int)((time / deltaT) + 1);
}
int getStepsX(double deltaX){
  return (int)(((B - A) / deltaX) + 1);
}

// =======================================
// Дано:
double psi(double x){ // функция, задающая распределение электронной плотности в нулевой момент времени
  //return ((x >= 0) && (x <= 0.5)) ? cos(2*M_PI*x)*sin(2*M_PI*x) : 0;
  return ((x >= 0) && (x <= 0.5)) ? sin(2*M_PI*x) : 0;
}
double fi_a(double t){ // левое граничное условие
  return 0;
}
double fi_b(double t){ // правое граничное условие
  return 0;
}
double f(double currU, double x, double t){ // функция, задающая внешнее воздействие
  return (((x >= 0.58) && (x <= 0.6))||((x >= 0.78) && (x <= 0.8))) ? 5 : 0; //2 : 0;
  //return ((x >= 0) && (x <= 2)) ? 2 : 0;
}

// =======================================

double **calculation( // расчет электронной плотности с помощью неявной разностной схемы
                      double deltaX, // шаг по X
		      double deltaTShow, // шаг во времени, с которым запоминаем
		      double time // кол-во моделируемого времени
		    ){
  double deltaT = deltaTShow;
  // считаем коэффициенты a, b, c
  // WARNING пердпологаем что коэффициенты не зависят от координат и времени
  // хотя в статье при коэффициентах были индексы по x
  double coffA = (-(SIGMA * deltaT) / (deltaX * deltaX));
  double coffB = (1 + 2 * ((SIGMA * deltaT) / (deltaX * deltaX)) + K * deltaT);
  double coffC = coffA;
  
  // определяем количество шагов
  int stepsX = getStepsX(deltaX);
  int stepsT = getStepsT(deltaT, time);
  
  // выделяем память для результата
  double **data = new double*[stepsT];
  for(int i = 0; i < stepsT; i++){
    data[i] = new double[stepsX];
  }
  // Применяем начальные условия
  data[0][0] = fi_a(0);
  data[0][stepsX - 1] = fi_b(0);
  for(int x = 1; x < stepsX - 2; x++){
    data[0][x] = psi(x * deltaX);
  }
  
  // два массива на смену для коэффициентов прогонки
  double **alpha = new double*[2];
  double **betta = new double*[2];
  for(int i = 0; i < 2; i++){
    // от 0 до N-2 ~ 1 до N-1
    alpha[i] = new double[stepsX - 1];
    betta[i] = new double[stepsX - 1];
  }
  
  // непосредственно вычисления	
  double ksi;
  
  for(int t = 1; t < stepsT; t++){
    // определяем начальные условия
    alpha[t % 2][0] = 0;
    betta[t % 2][0] = fi_a(t + 1);S
    
    data[t][0] = fi_a(t);
    data[t][stepsX - 1] = fi_b(t);
  
    // прямой ход метода прогонки
    for(int x = 1; x < stepsX; x++){
      alpha[t % 2][x] = coffA / (coffB + coffC * alpha[t % 2][x - 1]);
      ksi = data[(t - 1) % 2][x] + deltaT * f(data[(t - 1) % 2][x], x, t);
      betta[t % 2][x] = (ksi - coffC * betta[t % 2][x - 1]) / (coffB + coffC * alpha[t % 2][x - 1]); 
    }
    
    // обратный ход метода прогонки
    for(int x = stepsX - 1; x > 0; x--){
      // расчитываем непосредственно электронную плотность на t шаге по времени
      data[t][x] = alpha[t % 2][x] * data[t][x + 1] + betta[t % 2][x];
    }
  }
  
  
  return data;
}

int main(int argc, char **argv){ 
  double deltaX = atof(argv[1]); // шаг сетки по х // в статье h
  double deltaTShow = atof(argv[2]); // шаг по времени в результате
  double deltaT = deltaTShow;
  // окончание Show осталось исторически от явного метода =)
  double time = atof(argv[3]); // время моделирования
  
  
  double **data = calculation(deltaX, deltaTShow, time);
  
  
  int stepsX = getStepsX(deltaX);
  int stepsT = getStepsT(deltaT, time);
  // вывод в JSON:
  // {
  //   u: [[]] // функция распределения электронной плотности
  //   f: [[]] // функция внешнего воздействия\
  // }
  std::cout << "{\"u\":[";
  for(int i = 0; i < stepsT; i++){
    std::cout << "[";
    for(int j = 0; j < stepsX; j++){
      std::cout << data[i][j];
      if(j < (stepsX - 1)){
        std::cout << ",";
      }
    }
    std::cout << "]";
    if(i < (stepsT - 1)){
      std::cout << ",";
    }
  }
  std::cout << "],\"f\":[";
    for(int i = 0; i < stepsT; i++){
    std::cout << "[";
    for(int j = 0; j < stepsX; j++){
      std::cout << f(data[i][j], j * deltaX, i * deltaTShow);
      if(j < (stepsX - 1)){
        std::cout << ",";
      }
    }
    std::cout << "]";
    if(i < (stepsT - 1)){
      std::cout << ",";
    }
  }
  std::cout << "]}";
  return 0;
}